x_mesh:
[[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]]
y_mesh:
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[2. 2. 2. 2. 2. 2. 2. 2. 2. 2.]
[3. 3. 3. 3. 3. 3. 3. 3. 3. 3.]
[4. 4. 4. 4. 4. 4. 4. 4. 4. 4.]
[5. 5. 5. 5. 5. 5. 5. 5. 5. 5.]
[6. 6. 6. 6. 6. 6. 6. 6. 6. 6.]
[7. 7. 7. 7. 7. 7. 7. 7. 7. 7.]
[8. 8. 8. 8. 8. 8. 8. 8. 8. 8.]
[9. 9. 9. 9. 9. 9. 9. 9. 9. 9.]]
Tutorial 9: Advanced Modelling and Visualisation
Programming for Engineers (ENGG1001)
Tutorial Tasks
Task 1: Creating a Meshgrid for 3D Surface Plots
So far, you have had practice with creating 1D numpy arrays and plotting them. For 3D surface plots, we need to create a grid of points in the x-y plane. Write a function create_meshgrid(length: float, num_points: int) -> tuple[np.ndarray, np.ndarray] that creates a meshgrid of points from -length to length with num_points for both x and y axes. For a further explanation of what a meshgrid is, see the callout below.
What exactly is a meshgrid? A meshgrid is a 2D grid of points that can be used for plotting surfaces or evaluating functions over a 2D domain. It consists of two 2D arrays. Consider the example below where we create a 10 x 10 meshgrid:
Essentailly, a meshgrid is a way of presenting a set of cartesian coordinates in a structured format. So if we feed an index (which essentially acts as a coordinate) into the meshgrid, we can get the corresponding x and y values. For example, if we want to get the x and y values at the index (2, 3), we can do:
x value at x_mesh[2, 3]: 3.0
y value at y_mesh[2, 3]: 2.0
Here is a sample usage of the function:
What function could you use to create the 1D arrays for x and y? Once you have created the 1D arrays, then you can use numpy.meshgrid to create the 2D grid.
Task 2: Creating a Plane Surface
Now that we can create a meshgrid of points, we can use it to create a 3D surface. Write a function plane_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, a: float, b: float, c: float) -> np.ndarray that computes the z values for a plane surface defined by the equation: \[z = ax + by + c\]
Here is a sample usage of the function:
To compute the wave surface, you will need to calculate the radial distance R for each point in the meshgrid and then apply the RBF formula. Remember to use vectorisation!
Task 3: Modelling waves with a Radial Basis Function
Now that we have a meshgrid of points, we can model a wave surface using a Radial Basis Function (RBF). The RBF is defined as: \[f(x, y) = A \cdot \cos(kr - \omega t) \cdot e^{-\alpha r}\] Where:
- \(A\) is the amplitude of the wave
- \(k\) is the wave number
- \(r\) is the radial distance from the origin: \(r = \sqrt{x^2 + y^2}\)
- \(\omega\) is the angular frequency
- \(\alpha\) is the decay constant of the wave with respect to distance
- \(t\) is time in seconds
Write a function radial_wave(x_mesh: np.ndarray, y_mesh: np.ndarray, A: float, k: float, omega: float, alpha: float, t: float) -> np.ndarray that computes the wave surface values for the given parameters. Here is a sample usage of the function:
Task 4: Visualising 3D Surfaces
Now that we have a function to compute the wave surface, we can visualise it using Matplotlib’s 3D plotting capabilities. Write a function plot_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, surface: np.ndarray) that creates a 3D surface plot of the surface.
Here is a sample usage of the function potting the plane surface:
Which produces the following 3D surface plot:
Also try plotting the wave surface using the same function:
Which produces the following 3D surface plot:
To be able to setup an axis for a 3D plot, you need to specify the projection argument when creating the subplot. The plot_surface method can be used to create the surface plot. You can use the code below to get started:
If you want to make your plot look nicer, you can also specify a colormap using the cmap argument in plot_surface. For example, you can use cmap='viridis' to use the Viridis colormap, like the example above.
Now that you have created an easy way to visualise the wave surface, determine what effects the following parameters have on the wave surface by changing their values and observing the resulting plot:
- Amplitude
A - Wave number
k - Angular frequency
omega - Decay constant
alpha - Num points in the meshgrid
num_points
Task 5: identifying local maxima and minima
Now that we have a visualisation of the wave surface, we can identify the local maxima and minima of the surface. Write a function find_local_extrema(surface: np.ndarray) -> tuple[np.ndarray, np.ndarray] that finds and returns the radius values (how far from the origin) of the local maxima and minima of the surface. Here is a sample usage of the function:
To find the local maxima and minima of radial function, we can use the fact that the surface is “radially symmetric”, meaning that the surface values only depend on the radius from the origin. As such, we can compute the radius values for a central “slice” of the surface (e.g. a given row or column of the surface), and then find the local maxima and minima of that slice.
To illustrate this, consider the following central slice of the curve above:
We can see that the local maxima and minima of the surface correspond to the local maxima and minima of this slice. Therefore, we can find the local maxima and minima of this slice to determine the radius values of the local maxima and minima of the surface.
To find the local maxima and minima of a 1D array, you can start by using the np.diff function to compute the differences between consecutive elements. Then you can find where the sign of these differences changes (i.e. from positive to negative for maxima, and from negative to positive for minima). Then you can store the corresponding radius values for those indices and return them as the output of the function.
Task 6: Animating a 3D Wave Surface
For a challenge task, try to create an animation of the wave surface over time. However there is a slight twist: the equation for the wave surface now accounts for time decay as well, so the wave surface will not only change with time, but it will also decay over time. The new equation is: \[f(x, y, t) = A \cdot \cos(kr - \omega t) \cdot e^{-\alpha r} \cdot e^{-\beta t}\] Where \(\beta\) is the decay constant of the wave with respect to time. You can use Matplotlib’s animation capabilities to create the animation. Here is a sample usage of the function:
Use the following incomplete code as a starting point for your animation function. You will need to complete the update_image function and the animate_wave_surface function to create the animation.
from matplotlib import animation
def update_image(frame, x_mesh, y_mesh, A, k, omega, alpha, beta, ax):
"""Updates the wave surface for each frame of the animation.
Parameters
----------
frame : int
The current frame number.
x_mesh : np.ndarray
The meshgrid of x coordinates.
y_mesh : np.ndarray
The meshgrid of y coordinates.
A : float
The amplitude of the wave.
k : float
The wave number.
omega : float
The angular frequency.
alpha : float
The decay constant of the wave with respect to distance.
beta : float
The decay constant of the wave with respect to time.
ax : matplotlib.axes._subplots.Axes3DSubplot
The 3D axis object to update.
"""
wave_surface = # write your code here!
ax.clear()
ax.plot_surface(x_mesh, y_mesh, wave_surface)
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Wave Surface')
ax.set_title('Animated 3D Wave Surface Plot')
return [ax]
def animate_wave_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, A: float,
k: float, omega: float, alpha: float, beta: float,
delta_t: float, num_steps: int):
"""Creates an animation of the wave surface over time.
Parameters
----------
x_mesh : np.ndarray
The meshgrid of x coordinates.
y_mesh : np.ndarray
The meshgrid of y coordinates.
A : float
The amplitude of the wave.
k : float
The wave number.
omega : float
The angular frequency.
alpha : float
The decay constant of the wave with respect to distance.
beta : float
The decay constant of the wave with respect to time.
delta_t : float
The time step for each frame of the animation.
num_steps : int
The number of steps in the animation.
"""
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
#| decorate the plot with labels and title here!
anim = animation.FuncAnimation(fig, update_image, frames=np.arange(0, num_steps, delta_t),
fargs=(x_mesh, y_mesh, A, k, omega, alpha, beta, ax), interval=50)
plt.show()